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ABSTRACT 

We present a technique for estimating the mass in the outskirts of galaxy clusters 
where the usual assumption of dynamical equilibrium is not valid. The method assumes 
that clusters form through hierarchical clustering and requires only galaxy redshifts 
and positions on the sky. We apply the method to dissipationless cosmological TV- 
body simulations where galaxies form and evolve according to semi-analytic modelling. 
The method recovers the actual cluster mass profile within a factor of two to several 
megaparsecs from the cluster centre. This error originates from projection effects, 
sparse sampling, and contamination by foreground and background galaxies. In the 
absence of velocity biases, this method can provide an estimate of the mass-to-light 
ratio on scales ^ 1 — lO/i^^ Mpc where this quantity is still poorly known. 



Key words: dark matter 
miscellaneous. 



1 INTRODUCTION 

At scales ^ Ih^^ Mpqj, dynamics of galaxy groups and 
of the central region of clusters give mass-to-light ratios 
M/ L ~ 20 — 30% of that required to close the Universe (e.g. 
Carlberg et al. 1996; Ramella, Pisani, & Geller 1997). Ob- 
servations of light curves of Type la supernovae at redshifts 
z < 1 confirm that there might not be indeed enough mass to 
close the Universe (Riess et al. 1998; Garnavich et al. 1998; 
Perlmutter et al. 1998). On the other hand, from the anal- 
ysis of velocity fields on mildly non-linear scales (^ lQh~^ 
Mpc), we can measure the quantity j3 = f2o'^/b (e.g. Strauss 
& Willick 1995) where Slo is the mean mass density of the 
Universe and 6, assumed independent of scale, is the bias 
parameter, the ratio between the galaxy and the mass over- 
density fields. Recent estimates yield /3 ~ 0.5 — 1.1 (e.g. 
Dekel, Burstein & White 1997; Schmoldt et al. 1999). If 
galaxies cluster more than mass (6 > 1), this result implies 
large values of f2o- 

At intermediate scales, ~ 1 — lO/i"^ Mpc, estimating 
the mass-to- light ratio is particularly difficult. These scales 
correspond to the immediate surroundings of rich clusters, 
the largest non-linear systems forming at the present epoch. 
Matter is falling onto the cluster for the first time. Therefore, 
neither dynamical equilibrium nor linear theory are valid 
descriptions of the dynamics of this falling matter: at these 
scales, clusters are already in the non-linear regime although 
not yet virialized. 



* We use the Hubble constant Hq 
throughout. 
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We can derive constraints on Q.o from the dynamics of 
these infall regions on a statistical basis, by modelling the 
galaxy-cluster cross correlation function (e.g. Croft, Dalton 
& Efstathiou 1999). On a single cluster basis, however, us- 
ing the spherical infall model to extract information on Q.a 
(Regos & Geller 1989) seems impractical, because in hierar- 
chical clustering scenarios random motions are an important 
component of the velocity field (van Haarlem & van de Wey- 
gaert 1993). Nonetheless, the galaxy distribution in redshift 
space still contains information about the cluster mass (Di- 
aferio & Geller 1997, DG hereafter). 

DG suggest a method for measuring the mass profile 
at radii larger than the virial radius of galaxy clusters. DG 
show that, if the redshift space coordinates of the dark mat- 
ter particles were measurable, we could estimate the mass 
profile with a ~ 30% accuracy, independently of the mass 
and the dynamical state of the halo. 

Here, we present an operational method for extracting 
the mass profile from the redshift space distribution of galax- 
ies within real clusters. This mass estimator relies on kine- 
matic data alone. Thus, the mass estimate is independent 
of the relative distribution of mass and light, unless galaxies 
are not good tracers of the velocity field in the infall re- 
gion. Detailed modelling suggests that velocity bias is very 
weak on these scales (Kauffmann et al. 1999a; Diaferio et 
al. 1999). We apply the method to catalogues of galaxies 
formed and evolved using semi-analytic procedures within 
the dark matter halos of dissipationless A'^-body simulations. 
We show that it is indeed possible to measure the mass of 
clusters to large radii. Because a few hundred galaxy red- 
shifts are required to give reliable results, this measurement 
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is feasible only on rather massive clusters. This method en- 
ables us to estimate the value of the mass-to-light ratio on 
scales ~ 1 - 10h~^ Mpc. 

Sections ti and H review the DG interpretation of the 
galaxy density distribution in the redshift diagram of clus- 
ters and the assumptions of their mass estimation method. 
Section W outlines the operational method of extracting the 
mass profile from the redshift diagram, and Section g shows 
some applications to galaxy clusters extracted from two cos- 
mological models (Section pi). In Appendix M we also de- 
scribe a method of locating the centre of a galaxy cluster 
from a list of galaxies with redshift space coordinates. 




2 INTERPRETING THE REDSHIFT 
DIAGRAM 

Observations provide three out of the six phase space co- 
ordinates of a galaxy. Define the redshift diagram as the 
plane {r,v), where r is the galaxy angular separation from 
the cluster centre and v is its line-of-sight velocity relative 
to the cluster centre of mass. 

The spherical infall model predicts the existence of two 
curves in this plane where the galaxy number density is in- 
finite (Regos & Geller 1989). These caustics define a char- 
acteristic "trumpet" shape and enclose galaxies which are, 
in real space, at distances both smaller and larger than the 
turnaround radius from the cluster centre. Galaxies outside 
the caustics are only at distances larger than the turnaround 
radius. We define the amplitude A{r) of these caustics in 
redshift space as half the difference between the upper and 
the lower caustic. The sphericall infall model predicts a de- 
pendence of A{r) on the cosmological density parameter f2o 
(Regos & Geller 1989; DG). 

A clean contrast between the regions interior and ex- 
terior to the caustics seem indeed evident in the redshift 
diagrams of some real clusters (Quintana, Ramirez, & Way 
1996; den Hartog & Katgert 1996; Geller, Diaferio & Kurtz 
1999a; Geller et al. 1999b). However, it is a challenge to 
identify the location of the caustics reliably (e.g. Regos & 
Geller 1989; van Haarlem et al. 1993). Van Haarlem & van 
de Weygaert (1993) were the first to use A''-body simula- 
tions of dark matter halos formed in hierarchical clustering 
scenarios to show that (1) even in redshift diagrams of dark 
matter particles the caustics can be poorly determined, and 
(2) the amplitude of the plausible caustics is usually larger 
than the amplitude predicted by the spherical infall model 
for a given cosmology. They concluded correctly that the 
identification of caustics in cluster redshift diagrams is not 
a reliable method of estimating Qq ■ 

Van Haarlem & van de Weygaert (1993) did not clar- 
ify, however, what ultimately produces the caustics in the 
redshift diagram of real clusters. DG suggest a heuristic ar- 
gument to explain both the presence and the amplitude of 
the caustics. In hierarchical clustering scenarios non-radial 
motions and substructure obscure the pure radial motion 
predicted by the spherical infall model and infiate the pre- 
dicted amplitude of the caustics. How much do random mo- 
tions inflate this amplitude? In the infall region of a cluster, 
a galaxy with a large enough velocity can escape the cluster 
gravitational field in a time shorter than the age of the Uni- 



Figure 1. Geometry of the problem. An observer O views the 
cluster centre C along the line of sight OC and the galaxies A 
and B along the line of sight OAB at angular separation ip from 
the cluster centre. The longitudinal angles of A and B are 6 + -f 
and 8, respectively. 



verse.n We may therefore expect only a few escaping galaxies 
lying outside the caustic region in the redshift diagram; this 
fact rapidly decreases the galaxy number density outside the 
caustics. Thus, the appearence of the "caustics" is related 
to the local escape velocity rather than to the cosmological 
density parameter. In other words, the caustic amplitude is a 
measure of the cluster gravitational potential at a particular 
distance from the cluster centre. 

Of course, we do not have a correlation between caustic 
amplitude and cluster gravitational potential when a galaxy 
system falling onto the central cluster has a mass comparable 
to that of the cluster; in this case, the amplitude of the 
caustics is related to the gravitational potential of this falling 
system rather than the cluster itself. 

In order to relate the caustic amplitude A{r) to the clus- 
ter gravitational potential (f>{r) we proceed as follows (DG).n 
Regardless of the stability of the cluster, at any particular 
shell of radius r the maximum allowed velocity is the escape 
velocity Vg{r) — —2(j){r). The line-of-sight component of this 
escape velocity sets the location and the amplitude of the 
caustics in the redshift diagram. 

We now compute this expected maximum observable 
velocity. The line-of-sight component of Vs{r) at any given 
angular separation r± from the cluster centre is (see e.g. 
point B in Fig. |I|) 

T As a rough estimate of the astrophysical quantities involved, 
consider the gravitational field generated by a point of mass M. 
The escape velocity at a distance r from M is 

^o(r)-927 ,/" r kms-l. (1) 

In such a field, the deceleration is 



(2) 



v(r) ~ -440 I , ' , I ( ^ I h km s^^Gyr^^ 

Recall that a velocity of 1000 km s^^ corresponds to ^ 1 Mpc 

Gyr-i. 

T We restrict our argument to spherical symmetry because we 

expect the azimuthal average to enhance the appearence of the 

caustics in the redshift diagram. 
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«c,ios(r|rx) = «c,e(r) 



r± 



(r) 



ir' 



vl/2 



(3) 



where the longitudinal angle d is the angle between 
the galaxy position vector r centered on the cluster 
and the line-of-sight to the cluster centre. The equation 
dVe^ios{'''\r±)/dr = yields the maximum observable veloc- 
ity along the line-of-sight at fixed angular separation r±. 
The radius r = rM, where this maximum occurs, generally 
differs from rj_ and depends on the velocity profiles Ve,r{r) 
and Vc,e(r), i.e. the gravitational potential (j){r) and the tidal 
field determining the non-radial component of the velocity 
field. 

We seek an expression for the expected amplitude 
A{r±) of the caustics which is independent of the form of 
the escape velocity profiles. Consider the cluster velocity 
anisotropy parameter /3(r) = 1 — {vg+v^)/2{v'l:), where v^ is 
the azimuthal component of a galaxy velocity v, and the an- 
gle brackets indicate an average over the velocities of all the 
galaxies within the volume d"^r centered on position r. If the 
cluster rotation is negligible, {vg +v^) — 2{vj) = (n^) — {v^) 
and we can write {v'^{rj_)), the mean of the square of the 
escape velocity on the shell of radius r±, by combining its 
line of sight component v^i^^{r±) = v'^ g{r±) (equation yj) 
and /3{r±): 



{fc,9{r±)) = (we(rx) )- 



-/?(r-i 



3 - 2/3(rj 
Now, define the quantities 

3-2/3(rx) , . . 



9iP) 



{r± 



2l<^(rx)l. 



(4) 



(5) 



l-/?(rx)' ---'- g(p) ' 

g{P) contains the information about the anisotropy of the ve- 
locity field governed by the infall and the tidal field; (/f>/3(rx) 
combines g{f3) with the potential 4i{r±) determined by the 
mass distribution alone. By assuming 2|(/>(rx)| — {Vg{r±)), 
we can rewrite equation (vw as 



(j}f3{r±) ~ {vlg(r±)) = {We,los(r±)>. 



(6) 



When rM ^ r±, v^i^^ir = rM\r±) > v^i^^ir = rx|rx). 
Thus, a priori, the true amplitude (Wc,ios('' = rM\r±)) of 
the caustics is larger than 4'i3{r±). Moreover, the measured 
A{r±) can be larger or smaller than («o,ios('~ = rM\r±)) 
depending on the presence or absence of escaping galaxies 
in the redshift diagram. However, comparisons with clus- 
ters in Ai'-body simulations indicate that (jypir^), computed 
with the full phase space information, agrees, within the un- 
certainties, with the amplitude J\?(r±^ computed from the 
redshift diagram (DG; see also Section H). We can therefore 

fri) 



consider A? (rx ) as a measure of i 

Mri-)^A\r^). 



(7) 



The equation above represents our physical interpretation 
of the caustic amplitude A{r_i_) in the redshift diagram of 
clusters. 

Equation (H) clearly shows that non-radial motion is the 
fundamental ingredient of this interpretation of the caustic 
amplitude A{r±). In fact, if orbits were radial, /3(rx) = 1, 
and equation (Bl) would predict an uncorrect (^^(rx) = 0. 
On the other hand, random motions dominate the velocity 
field of the virialized region of the cluster. Thus, this inter- 
pretation holds both in the central and in the infall region 
of the cluster. 



In conclusion, in redshift diagrams of real systems, we 
expect a concentration of galaxies around the cluster red- 
shift at any fixed angular separation. Because of infall, we 
expect a somewhat clear separation between falling galax- 
ies and unrelated galaxies which are more distant than the 
turnaround radius. Because of random motions, this sepa- 
ration is not as sharp as predicted by the spherical infall 
model, but it is still apparent, and it is indeed observed 
(Geller et al. 1999a). Finally, the cluster gravitational po- 
tential determines the amplitude of the caustics. 



3 ESTIMATING THE INTERIOR MASS 

Suppose we have a method to measure A{r) and g{P). By 
assuming A?{r) = (jjpir), equation (0) readily yields, for the 
cumulative cluster mass GM{< r) — r'^d(j>/dr, 



GM{<r) = --A'{r)g{l3) 



d\nA^ ding 



dlnr 



+ 



dlnr 



(8) 



The two logarithmic derivatives in equation (fel) are compa- 
rable. Recently, Merrifield (1998) has shown that P{r) could 
actually be measured by detecting N ^ 50 galaxy wakes in 
the X-ray emission of clusters. Bent lobes of radio galaxies 
might also be used to cross-check these measurements; these 
galaxies will not be useful for actual measurements, however, 
because the expected number of these galaxies is not larger 
than ~ 5 per cluster (Cress 1999) . In any case, even if we 
had a measure of l3{r), a serious problem occurs indepen- 
dently of the method we use to estimate A^{r): sparse sam- 
pling and background and foreground galaxies will yield a 
very noisy A^{r); therefore its differentiation is not practical 
even after substantial smoothing (DG). 

Nevertheless, a measure of A{r) alone can still be used 
to estimate the cluster mass. DG suggest casting the relation 
between A{r) and M{< r) in the form 



•Jrn 



GM{< r) - GM{< ro) = / A^{x)Tf3ix)dx (9) 

where 

,p(r)r2 3-2/3(r) 



Mr) 



~2ttG- 



(r) l-/3{r) 



nr)g{P), 



p{r) is the cluster mass density profile and 



GM{< r) 



-4nG 



p{x)xdx 



(10) 
(11) 

(12) 



is the gravitational potential generated by the cluster. 

Equation (H) integrates the estimated A? (r) rather than 
differentiating it, thus averaging out its un-physical rapid 
variations. However, the crucial property of equation (Bl) is 
that the function T/jir) is slowly varying at large radii in 
hierarchical clustering scenarios, as we will show below. This 
property is of major importance, because we may assume 
jr^(r) — const = T(i. Equation (H) then becomes 

IT 

GM{<r)~GM{<ro)=Ti3 A\x)dx. (13) 

Thus, we make the mass estimation method independent of 
any scale length, implicitly assumed in the function Ti3{r-), 
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Figure 2. Filling function .^nkwC*") for the NFW density pro- 
file (equation h4). The bold line is .7^nfw(''/^s)- Thin lines are 
■^NFw('') for different values of the concentration parameter 
c = r-20o/''s- 



although still dependent on the unknown . 7-"^. Of the three 
equations (pi), (bl) and (|l3), only equation (hsl) turns out to 
be of practical relevance. We emphasize that equation ( |l3[ ) 
is justified only in hierarchical clustering scenarios. 

Note that the similarity between equation y^ and equa- 
might suggest using an iterative procedure: equa- 
provides the first mass profile from which we can 
estimate Tf3{r) and we can then use equation (tel) until con- 
vergence is achieved. Unfortunately, J-/3{r) contains infor- 
mation on both the mass distribution and the velocity field 
(equation hll), whereas equations (bl) and (113) return infor- 
mation on the mass profile alone. Therefore, we should still 
make an assumption about g{P)- However, even in this case, 
the improvement in the mass estimate is likely to be insuffi- 
cient to reduce the uncertainty which is introduced by pro- 
jection effects. These effects actually dominate the problem, 
as we will show in Section pi 

To justify J-/3{r) = T(r)g{l3) — const in hierarchical 
clustering cosmogonies, first consider the function J-{r) = 
—2-KGp{r)r^/(f){r), and a density profile approximated piece- 
wise by power laws p{r) oc r^°'^''' . Provided the mass is 
finite at r = {a{r) < 3, when r -^ 0), and (/)(r) con- 
verges at large radii {a{r) > 2, when r -^ oo), we obtain 
T{r) oc r2-"('-)/[M(< r)r-'^ + r^-°'^% At large radu, T{r) 
behaves differently depending on the behavior of Af (< r) oc 
^3-a(r) j^ ^j^g rnass is finite (a(r) > 3, when r -^ oo), T{r) 
decreases as a power law, T{r) oc r^'"^^'^'. If the mass di- 
verges (2 < a{r) < 3, r ^ oo), J-{r) ~ const. This is the 
case of interest here. 

A'^-body simulations of hierarchical clustering scenar- 
ios suggest a universal density profile for dark matter halos 
(Navarro, Frenk & White 1997, NFW hereafter). This pro- 
file yields a diverging mass and a convergent gravitational 
potential at the same time, as required by the argument 
above. However, in this case, a{r) = 3 when r —> oo, imply- 
ing M(< r) oc ln(r) and 



where r^ is a scale length defined in terms of the halo con- 
centration c — r2oo/f's, and r2oo is the radius of the sphere 
whose average mass density is 200 times the critical density. 
J'NFwij') is not constant in hierarchical clustering scenarios 
(see Figure 0) . However, .Tt^fw if) decreases by ~ 20% at 
most in the radius range r G [1, 10]r2oo for c £ [1, 7]. This 
range of concentration c is typical for massive halos in Cold 
Dark Matter (CDM) models; in fact, in our A-body mod- 
els (Section H), 90% of halos more massive than lQ^'^h~^MQ 
have c in the range [1.6, 5.2] with median c = 3.7. 

To see that J-0{r) = J-{r)g{P) is also a slowly vary- 
ing function of r, we compute the functions /3(r) and g(d) 
for massive clusters in our A-body simulations (Section m. 
Because these two functions do not show strong variations 
with r, the assumption of a slowly varying Tfj{r) appears 
reasonable. We will specify the value of the constant Tp in 
Section ra. 

Finally, it is important to note that the identification 
of the amplitude A{r) of the caustics with the gravitational 
potential [4'I3{''')Y at any radius r is independent of the den- 
sity profile and the dynamical state of the system, provided 
that random motions contribute significantly to the velocity 
field and the spherical assumption is a good approximation. 
On the other hand, the measurement of Air) yields the mass 
profile M(< r) (equation |l3| ) only when the density profile 
is proportional to r""'"^', with 2 < air) < 3 when r — > oo. 



4 MEASURING A{r) 

We now outline an operational method for locating the caus- 
tics in the redshift diagram of real clusters. Compiling a red- 
shift diagram requires knowledge of the cluster centre and 
its radial velocity. Appendix vM describes how to determine 
the cluster centre when only the celestial coordinates of the 
galaxies and their redshifts are available. 

In order to determine the amplitude A{r) of the caustics 
at fixed r, consider the two-dimensional density distribution 
function fq{r,v)drdv, namely the number of galaxies with 
projected separation in the interval (r, r + dr) and line-of- 
sight velocity in the interval {v, v+dv). In the spherical infall 
model, at fixed r, fq{r,v) reaches infinity at the two caustic 
locations. Random motions, however, wash these two spikes 
out and fq{r, v) has a maximum close to the cluster redshift; 
in the absence of massive substructure and foreground and 
background galaxies, and in the absence of escaping galaxies, 
the solutions of the equation fq{r,v) = would determine 
the amplitude A{r). However, this situation never occurs 
because, even for isolated systems, escaping galaxies can al- 



ways be present. Thus, fq 



always becomes zero outside 



.7T^Fw(r) — 



1 jr/rs)^ 1 

2(l + r/r,)2 1n(l + r/r,) 



(14) 



the actual location of the caustics. 

We therefore need a recipe for choosing a threshold k 
such that the equation fq{r,v) = k determines the ampli- 
tude A{r). Moreover, we have to face the non-trivial issue 
of accurately estimating /^(r, u) where it is close to zero. It 
is clear that the this task is not easy, because sparse sam- 
pling leads to an underestimate of fq{r,v) and the presence 
of foreground and background galaxies leads to an overes- 
timate of fq{r,v) at the caustic loc atio n. We describe the 
determination of fq{r,v) in Section 4.1, In Section 4.2 we 
describe the choice of k. 
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4.1 Estimating the redshift diagram density 
distribution 

Consider A*' galaxies with coordinates x = (r, v), where r and 
V are conveniently rescaled (we come to this issue later). We 
use an adaptive kernel method (Silverman 1986; see also 
Pisani 1993; 1996) to estimate the density distribution of 
galaxies within the redshift diagram 



/. 



M^^i:i^-(^) 



where 



Kit)^{l--'(^ 



otherwise 



(15) 



(16) 



and hi = hchoptXi is a local smoothing length depending on 
the local density. The optimal smoothing length is 



/^op 



3.12 f (j^ -\~ (jy 



iVl/6 



(17) 



where a,, and a^ are the marginal standard deviations of 
the galaxy coordinates. The local smoothing factor is A; — 
['y/fi-{^'-)Y where /i is equation ( [L5[ ) where he — Xi = 1 
for any i, and log 7 = J2i ^oglfi (x,)]/A''. 

The average degree of smoothing is controlled by he- 
The optimal he should minimize the integrated square error 
between the estimator /q(x) and the true (unknown) density 
/(x), e{hc) = J[fq{x) — /(x)]^d^x. However, it is easy to 
show that minimizing e{hc) is equivalent to minimizing 



Mo (ft.) = 



N 



f^(^)d^^-^y-f^~ 



2 



'(x.) 



(18) 



f(-')/ 



where fq (xi) is the density estimated at Xi using all the 
data points except Xi (Silverman 1986). Therefore, Moihc) 
can be estimated with the data alone and we do not need to 
assume any form for the true density /(x). 

We now address the question of rescaling r and v, such 
that we can use spherical smoothing windows of size hi. 
The adaptive kernel method is designed to estimate density 
distributions of random variables. If this were the case, we 
could apply the whitening transformation, namely we would 
linearly transform the data to have a unit covariance matrix 
(Fukunaga 1990). Unfortunately, x = (r, v) is not a random 
variable and the elements of the covariance matrix depend 
on the limits imposed a priori on the redshift diagram. 

We therefore have to rescale r and w in a sensible way. 
Let us use the Hubble constant to have r and v in the same 
units. It is obvious that, for example, 100 km s~^ in the v 
direction shold not have the same weight as lh~^ Mpc = 100 
km s~^ in the r direction. We thus rescale r and v such 
that the ratio h^/hr of the smoothing window sizes along v 
and r respectively takes a chosen value q. Galaxy redshifts 
have typical uncertainties of 50 km s~^. Positions of galaxies 
within nearby clusters have uncertainties of ~ Q.Q2h^^ Mpc. 
We thus set q — 25 and keep this value fixed hereafter. Note 
however that difi'erent values of q in the range [10, 50] have 
little effects on the results. 



4.2 Choosing the threshold 

We determine the amplitude A{r) at fixed r by finding 
the solutions of the equation fq{r,v) = k. Specifically, the 
first upper and lower solutions Vu and Vd, away from the 
maximum of fq (r, v) closest to v = 0, determine the am- 
plitude A{r) = mindful, |vd|}. Note that the prescription 
A{r) = Tnm{\vu\, \vd\} is equivalent to A{r) = {vu — «d)/2 
for an isolated spherically symmetric system. However, our 
prescription is more robust than A{r) = {vu — Vd)/2 against 
interloper contamination and the presence of massive sub- 
structure. 

It is clear that, although we have determined fq{r,v) 
uniquely (except for the choice of q), there are an infi- 
nite number of thresholds k we can use to determine A{r). 
It seems reasonable to assume that, at least in the cen- 
tral region, the cluster has attained virial stability. There- 
fore, in this region, the equation (Wesc)-R = 4(ii^)_r must 
hold, where now the angular brackets indicate an average 
over the whole sphere of radius R, and velocities are three- 
dimensional. In the data we only have one-dimensional in- 
formation, and we thus further assume that, if the veloc- 
ity field is roughly isotropic in the central region, our ex- 
pression also holds when v is the galaxy line-of-sight ve- 
locity and (i'csc)k,-R ~ L -^^{''')vi''')d-r/ J, (p{r)dr, where 
(p{r) = J fq{r,v)dv; {v^sc)k,,R is the only K-dependent quan- 
tity. We can choose re by minimizing the function 

S{K,R) = \{vL),.,R-A(v^)Rf. (19) 

In Appendix |A| we describe a procedure for locating the clus- 
ter centre from a list of galaxy positions. This procedure also 
identifies the cluster members. In the following, we define R 
as the mean projected distance of the members from the 
cluster centre, and {v'^)r the one- dimensional velocity dis- 
persion of the cluster members. 

The amplitude A{r) is obviously sensitive to the value 
of K. If fq{r,v) has a rather shallow gradient towards the 
minimum around the real location of the caustics, a slightly 
incorrect estimate of {v'^)r can lead to particularly inac- 
curate estimates of A{r). Thus, in these cases, minimizing 
S{n, R) will only give a starting value for the choice of k. The 
final value of k can be chosen by hand and the location of 
the caustics will unfortunately be subjective. This freedom 
is unsatisfactory. On the other hand, k remains a useful one- 
parameter tool which quantifies the complexity of the infall 
region of individual clusters and the contamination of their 
redshift diagrams by background and foreground galaxies. 
In any case, this subjective tuning turns out to be necessary 
in a few cases only, and mainly for one of the two cosmolog- 
ical models we investigate here, which yields the poorest fit 
to observation. This result suggests that the minimization 
of S(k,R) is likely to be a robust procedure for real clusters, 
as the Coma cluster has indeed already shown (Geller et al. 
1999a). 

To control the contamination by background and fore- 
ground galaxies efficiently, we need a final step. Consider 
the logarithmic derivative 2d\nA/dlnr — dln^^/dlnr; we 
have dliKpi^/dlnr = din |</f>j/dlnr — dlng/dlnr. For any re- 
alistic system dln\(f)\/dlnr ^ 0. Our A^-body simulations 
show that dlng/d\nr ^ —1/2. Thus, we can safely claim 
that din (p/}/ din r JS 1/2 should hold for any r. To be con- 
servative, we therefore accept only values of A{r) which yield 
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Parameters of the two GIF simulations (Kauffmann et al. 1999a) 
used in this paper. The Hubble constant h, the particle mass rrip, 
and the comoving size L of the simulation box are in units of Hq = 
100 km s~i Mpc^i, 10^" h^'^ Mq, and h~'^ Mpc, respectively. 

d\nA/dlnr < 1; otherwise we impose a new value of A{r) 
which yields 2d\nA/d\nr = 1/2. 



5 NUMERICAL MODELS 

In Section o we apply our mass estimation method to dis- 
sipationless A''-body simulations where we form and evolve 
galaxies within their dark matter halos with phenomeno- 
logical recipes. Here, we briefly describe our cosmological 
simulations. 



5.1 The GIF simulations 

We use two TV-body simulations from the GIF project 
(Kauffmann et al. 1999a). These simulations use 256^ dark 
matter particles to model the evolution, from redshift z = 50 
to the present, of the dark matter density perturbations of 
a CDM universe with an initial power spectrum 

\2il.l3i -2/1.13 



P{k) oc fc{l+[6.4fc//ir+(3fc//ir)^''^+(1.7fc//ir)^]^-^^}" 



(20) 



where k is in units of h Mpc~^, and F — 0.21 is a shape 
parameter (Efstathiou, Bond & White 1992). The normal- 
ization is fixed by the ratio af of the variances of the mass 
and galaxy fiuctuations within randomly placed spheres of 
radius 8h~^ Mpc. The GIF models are normalized to give the 
correct abundance of rich galaxy clusters at the present time. 
The simulations were run with Hydra (Pearce & Couch- 
man 1997), the parallel version of the AP^M code (Couch- 
man 1991; Couchman, Thomas, & Pearce 1995), kindly 
provided by the Virgo supercomputing consortium (Jenkins 
et al. 1997). The two models we consider here are a flat 
model with or without a cosmological constant: Qq = 0.3 
and i^A = 0.7 (ACDM) and ^o = 1 (rCDM). Table § sum- 
marizes the parameters of these models. 

Kauffmann et al. (1999a) combine these A''-body sim- 
ulations with semi-analytic modelling to form and evolve 
galaxies within dark matter halos. They provide a detailed 
description of this procedure and compare the observable 
properties of the simulated galaxy catalogues with the real 
Universe in a series of papers (Kauffmann et al. 1999a; 
Kauffmann et al. 1999b; Diaferio et al. 1999) . The relevant 
physical processes for galaxy formation include gas cool- 
ing, star formation, supernova feedback, stellar evolution, 
and merging of galaxies. Previous attempts to compare such 
dissipationless simulations with the real galaxy distribution 
were based on some high-peak statistical model for galaxy 
formation or on ad hoc assumptions about the mass-to-light 
ratio of dark matter halos (e.g. Frederic 1995; Nolthenius, 
Klypin & Primack 1997). Alternatively, we could consider 
a full A'^-body/hydrodynamic simulation. However, the cur- 
rent state-of-the-art for such simulations can model only a 



small volume of the Universe, otherwise their resolution limit 
is larger than the size of galaxies and their halos (e.g. Frenk 
et al. 1996; Navarro & Steinmetz 1997; Weinberg, Hern- 
quist & Katz 1997; Gen & Ostriker 1999; Blanton et al. 
1999; Pearce et al. 1999). 

The GIF simulations have a resolution limit of 30/i~^ 
kpc and are suitable for investigating the reliability of our 
mass estimation method. We can in fact construct mock 
redshift diagrams with galaxies, along with their observable 
properties, rather than with dark matter particles alone. 
Note that, in these models, galaxy samples show only a weak 
velocity bias, independently of galaxy luminosities, as ex- 
pected if gravity is the driving force on scales ^ lh~^ Mpc. 

As discussed in Kauftmann et al. (1999a), the free pa- 
rameters entering the galaxy formation recipe have substan- 
tial effects on galaxy properties. Here, for consistency, we 
consider the "fiducial" models which yield reasonable fits 
to many but not all the observed properties of galaxies; in 
particular the model luminosity function is a rather poor fit 
to observation. For the purpose of this paper, it is worth 
noting that the rCDM model has a luminosity density in 
the B-band which is a factor ~ 2 larger than in real galaxy 
samples, whereas the ACDM model has a luminosity den- 
sity a factor ~ 2 smaller. In order to quantify the effect 
of the luminosity function on the derived properties of the 
galaxy distribution, Diaferio et al. (1999) assign new lumi- 
nosities to the model galaxies, while preserving their lumi- 
nosity rank in such a way as to reproduce the luminosity 
function of the Center for Astrophysics (CfA) redshift sur- 
vey (Marzke, Huchra, & Geller 1994) exactly. We therefore 
have two distinct sets of simulated galaxy catalogues: (1) 
the SALF catalogues where galaxies have luminosities de- 
rived from the semi-analyitic modelling, and (2) the CfALF 
catalogues where galaxies have luminosities imposed accord- 
ing to the GfA luminosity function. 

5.2 Profiles of dark matter halos 

Before testing the mass estimation method, we address the 
robustness of the assumptions described in Section ra, namely 
whether J-pir) is a slowly varying function of r in these CDM 
models. 

We identify dark matter halos using a friends-of-friends 
group finder which links particles closer than 0.2 times the 
interparticle separation. We take the position of the most 
bound particle as the centre of the halo. Here, we consider 
only dark matter halos with M200 > lO^^/i^^M©, where 
Af200 is the mass within r2oo- Less massive halos are less 
interesting for our purpose, because our mass estimation 
method, which requires the redshifts of a few hundred galax- 
ies, will be applicable to massive clusters only. 

Figure H shows the median profiles of the relevant quan- 
tities of the halos in both models. Shaded areas indicate the 
interquartile range of the profiles. The velocity anisotropy 
parameter /3(r) increases from ~ 0.2 to ~ 0.5 — 0.6 be- 
tween r — 0.1r2oo and r ~ 3r2oo, as radial motions be- 
come predominant. Remarkably, this result agrees with the 
value (3 ~ 0.5 estimated by Biviano et al. (1997) for the 
emission-line galaxies in the ESQ Nearby Abell Cluster Sur- 
vey. At the turnaround radius r « 4 — 5r2oo, (Vr) decreases 
faster than its non-radial counterparts (v^) and {vg); thus, 
/3(r) drops to negative values. Note that /3(r) is sensibly 
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Figure 3. Profiles of /3(r), g(/3), J^{r), and J^0ir) for the clusters with M200 > W^'^h-'^MQ in the rCDM (dashed lines) and ACDM 
(solid lines) models. Bold lines are the median profiles. Shaded areas show the interquartile range of the profiles. 



smaller than one everywhere, indicating that random mo- 
tions always contribute significantly to the velocity field, as 
expected in hierarchical clustering scenarios. The function 
(?(/3) we introduced in Section ti behaves as P{r). More im- 
portantly, g{f3) varies with r slowly, increasing by 30-50% at 
most between r = 0.1r2oo and r ^ 2 — 3r2oo- At radii where 
the Hubble flow dominates, /3(r) ^ 1 and g{P) -^ (x>. 

To compute the function J-{r) introduced in Section pi 
we need to compute the gravitational potential (ji{r) gen- 
erated by the halo. We could compute (ji{r) with equation 
([12) where we replace the upper limit of integration with 
a maximum radius. However, because the halo density pro- 
files are well approximated by the NFW formula, it is more 
efficient to obtain a fit to the density profile and compute 
(j}{r) analytically. This procedure bypasses the problem of 
subtracting the background at large radii. Figure Ja shows 
that J-{r) is slowly varying, as expected. Note that at radii 
»■ <; 3 — 4r2oo, the derivative of ^{r) is positive. In fact, the 
density profiles do not decrease as r~^ as predicted by the 
NFW profile, but rather decrease as r~°', with a ~ 1 — 2, 
because they must approach the background density. There- 
fore J-{r) ^ ~pr^/(f> increases as r"^~°/lnr. Note that the 
statistical spread in the y-{r) profiles accounts for different 
concentration parameters c. As expected, J^p (r) has a trend 
similar to the trend of J- if), because g(fi) does not vary 
strongly. Variations of (/(/3) only make the statistical spread 
of J'p{r) larger. 

Note that the difference between the rCDM and the 
ACDM models for each of the four profiles is mainly due to 
the different physical location of r2oo in the clusters of the 
two models: r2oo identifies a sphere with an average density, 
in units of the background density, which is larger in the low 
density universe than in the high-density one. 



In conclusion, ^pir) is not constant in these models, as 
assumed in equation (113) , but it varies by a factor of two at 
most within r ~ 3 — 4r2oo- Despite these variations, equation 
([13|) remains a useful recipe, as we show in the next section. 



6 GALAXY CLUSTERS 

We now apply our caustic location method and mass es- 
timation procedure to simulated clusters. To construct the 
redshift diagram, we use the centre of the dark matter halo 
and include galaxies brighter than Mb ~ — 18.5+51og/i. This 
magnitude limit corresponds to an apparent magnitude of 
m-iim = 15.7 at the distance of Coma {R = 7000 km s~^) and 
is one magnitude brighter than the limit of completeness of 
our galaxy catalogues (Kauffmann et al. 1999a). 

We also simulate observations of nearby clusters, by 
randomly choosing the observer's location on a sphere of 
radius R = 7000 km s~^ centered on the centre of the dark 
matter halo. We then compile a list of galaxies brighter than 
Tiiim = 15.7 with coordinates (a, 5, cz) within 10° from the 
halo centre; from this list, we locate the cluster centre with 
the method described in Appendix Ul This method usually 
locates the cluster centre accurately; both this procedure 
and the one which uses the centre of the dark matter halo 
directly, yield similar results. Below, we show results where 
we use the centre of the dark matter halo. 



6.1 Caustic Location 

Fig. y shows the redshift diagrams of a typical cluster in the 
ACDM model. The nine panels correspond to nine differ- 
ent lines of sight. The galaxies tend to populate a defined 
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Figure 4. Redshift diagrams of a cluster in the ACDM model viewed along nine different lines of sight. The circles show galaxies 
brighter than Mb = —18.5 + 51og/i. The bold lines show the identified caustic locations. When the contamination by foreground and/or 
background galaxies is too heavy for a reliable location of the caustics the method does not proceed further. This situation is represented 
by merging the caustics into a single horizontal line at the cluster redshift (see right-bottom panel). 



region in the redshift diagram. The contrast is quite evi- 
dent between the galaxies within and outside this region. In 
most cases, our method seems to locate the caustics, i.e. the 
borders of this region, correctly. 

As discussed in Diaferio et al. (1999), despite the fact 
that mock redshift surveys extracted from both the rCDM 
and the ACDM models do not show structures as sharply 
defined as in real surveys, the two models still show substan- 
tial differences: the ACDM model shows voids and filaments 
larger than the rCDM model and yields a better qualitative 
fit to the real Universe (see Schmalzing & Diaferio 1999 for 
a quantitative analysis). These differences also are quite ap- 
parent in redshift diagrams of clusters (Figs. Wand pi), where 
the density contrast between interlopers and galaxies falling 
onto the cluster is less evident in the rCDM than in the 
ACDM model. Nevertheless, in rCDM clusters, our method 
still locates the caustics reasonably well in most cases (Fig. 

The different redshift diagrams yielded by the two cos- 
mological models are related to the underlying mass distri- 
bution rather than to the galaxy luminosity function. Using 
the CfALF catalogue decreases the number of galaxies in the 
redshift diagrams of rCDM clusters but does not improve 
the appearence of the caustics. In fact, in this high density 
universe, clusters are still accreting mass at a large rate at 
redshift z — 0. On the other hand, in the ACDM, where 
using the CfALF catalogue increases the number of galaxies 
in the redshift diagrams, the accretion rate is substantially 



smaller and the caustics tend to appear more clearly, as ob- 
served in real clusters (Geller et al. 1999a). 

In principle, this redshift diagram morphology could 
be used to distinguish between high and low density uni- 
verses. For example, we could quantify this difference by 
using the distribution of the ratio between the number of 
galaxies lying within the caustics and the total number of 
galaxies within the redshift diagram. This ratio, of course, 
depends on the cluster and on the line of sight, but the 
moments of the distribution should be sensitive to the cos- 
mological model. For example, these ratios for the different 
lines of sight of the ACDM cluster (Fig. kh are in the range 
[0.32, 0.88] with median 0.53 (panel 2) and are considerably 
larger than the corresponding ratios of the rCDM cluster 
(Fig. |): [0.11,0.55] with median 0.39 (panel 7). These dis- 
tributions appear to be robust against variations of the mag- 
nitude limit of the galaxy sample. The feasibility of this kind 
of statistical tests will be investigated in future work. 

In panels 2, 5, and 9 of Fig. B the caustics seem to be 
outside the location we would guess subjectively, because 
of the heavy contamination by foreground and background 
galaxies. Our intuition is confirmed by comparing the ampli- 
tude A{r) of the caustics with the function [4'i3{''')Y com- 
puted from the full six-dimensional phase space information. 
The upper right panel of Fig. m shows all nine estimated pro- 
files for the rCDM cluster: the ones which tend to overesti- 
mate A{r) come from the three errant panels. In the other 
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Figure 5. Same as Fig. U for a cluster in the rCDM model. 
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Figure 6. Profiles of the two clusters shown in Figs. W and H. The thin lines are the profiles measured from the redshift diagrams. The 
bold lines in the four upper panels are the true profiles. 
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six redshift diagrams and in all the diagrams of the ACDM 
cluster the agreement is rather good. 

It is worthwhile to note, at this point, that our caustic 
location procedure is essentially an interloper-removal pro- 
cedure which uses the combined information on the galaxy 
position and line-of-sight velocity, unlike standard 3-sigma 
clipping procedures which ignore the position of galaxies 
(e.g. Yahil & Vidal 1977; Perea, del Olmo & Moles 1990). 
Our procedure enables us to estimate the mass profile of 
clusters as a by-product, on the assumption that clusters 
form through hierarchical clustering.^ 

Let us finally emphasize that our procedure does not 
remove all the interlopers: in fact, any interloper-removal 
procedure, which does not separate the galaxy redshift in its 
Hubble flow and peculiar velocity components, is obviously 
unable to identify interlopers lying within the caustics in 
redshift space. 



6.2 Mass Profiles 

To compute the mass profile from the caustic amplitude, we 
set Ti3 — 1/2 in equation dlq), as first suggested by DG, 
although Fig. H shows that the mean value of Tp may be 
slightly larger; moreover, we set ro = 0. We thus have 



GA/c 



t«0 = 2 



A {x)dx. 



(21) 



In writing equation (hsl) we have assumed that Tp{r) is 
roughly constant at large r. Thus, setting ro = is justi- 
fied if J-i3{r) is also roughly constant at small r and takes 
the same values as at a large r. In fact, Fig. 3 shows that 
J'fjir) £ [0.5, 1.0] for any r < 4r2oo, and Fig. 6 shows that 
the estimate of M(< r) at small r is as good as at large r. 

The mass profile is recovered within an uncertainty of 
50% out to ~ 8/i"^ Mpc for the ACDM cluster (left middle 
panel of Fig. |6|). For the rCDM cluster the estimated mass 
profile is less accurate at those distance, but still within a 
factor of two. 

Finally, the bottom panels of Fig. ra show the profiles of 
the _B-band mass-to-light ratio in units of the mean ratio for 
the entire simulation box. In the ACDM cluster, the mass-to- 
light ratio tends to be larger than the universal value at very 
small radii because of the deficiency of blue galaxies in the 
central region of clusters in this cosmology (Diaferio et al. 
1999). At radii ~ 4 — 6/i~^ Mpc both models show a mass-to- 
light ratio consistent with the universal value. Note that for 
these massive clusters, on average, the mass-to-light ratio 
within r2oo ~ 1.5/i~^ Mpc, computed with the full three- 
dimensional information, is close to the universal value (see 
Fig. 15 of Kauffmann et al. 1999a). This result is consistent 
with our Fig. ra where the estimation of the mass-to-light 



8 An interloper-removal procedure, which also uses both the po- 
sition and the velocity of galaxies, was introduced by den Hartog 
& Katgert (1996). Their iterative procedure is effective in deter- 
mining the cluster membership, but it cannot be used to estimate 
the cluster mass. In fact, this procedure estimates the profile of 
the maximum line-of-sight velocity allowed for a cluster member, 
namely a "caustic" , only after estimating the mass profile. More- 
over, the estimate of the mass profile relies on the virial theorem 
all the way out to the infall region. 



ratio at r ~ 1.5^~^ Mpc suffers from projection effects and 
errors in the mass estimate. 

Real clusters, of course, provide a single redshift dia- 
gram. The error in the measured value of A(r) should de- 
pend on the number of galaxies which contribute to the de- 
termination of A{r). We thus assume that the relative error 
5 Air) /Air) = K/max{/g(r, I))}, where the maximum value 
is found along the ii-axis at fixed r. We then define the er- 
ror in the mass profile as 8Mi = X/ -i |2™'j'5-4(rj)/^(rj)|, 
where rrij is the mass of the shell [rj-i, r_,] given by equation 
( [l3| ) with Tp = 1/2. This recipe yields errors in agreement 
with the typical spread due to the projection effects shown 
in Fig. 0. 



7 DISCUSSION 

DG suggested the possibility of using redshift data alone 
to measure the mass of clusters within ;$ 10/i~^ Mpc from 
their centres. This mass, combined with photometric mea- 
surements, provides an estimate of the mass-to-light ratio 
and therefore the mean mass density of the Universe, ilo, 
on the assumption that the value obtained on this relatively 
large scale is close to the global value. 

Here, we describe an operational procedure which can 
be applied to redshift diagrams of real clusters containing a 
few hundred galaxies with measured redshifts. We apply this 
procedure to galaxy clusters simulated from A''-body models 
which include semi-analytic modelling of galaxy formation. 
We are thus able to mock observations of clusters where 
the luminosity and the formation history of galaxies are in- 
cluded. We recover the actual cluster mass profile within a 
factor of two to several megaparsecs from the cluster center. 

For the sake of clarity we summarize here the assump- 
tions and the parameters entering our mass estimation pro- 
cedure, giving in parentheses the section where we discuss 
the issue extensively. We assume that (1) clusters are spher- 
ically symmetric (Section Q); (2) clusters form through hi- 
erarchical clustering: this assumption implies that (i) non- 
radial motions are an important component of the velocity 
field in the infall region of clusters (Section H), and (ii) the 
filling function Tf^ir), which combines the cluster density 
profile and the anisotropy of the velocity field, is roughly 
constant at large distance r from the cluster centre (Section 
bl); (3) substructures in the cluster infall region have mass 
substantially smaller than the cluster mass (Section bh; in 
other words, we cannot apply our mass estimation method 
when a major merging between clusters is taking place. 

The parameters governing the location of the caustics 
and the estimate of the mass profile are: (1) the ratio q be- 
tween the sizes of the smoothing window along the line-of- 
sight velocity v and the angular separation r (Section |4.l[ ); 
(2) the threshold k fo r th e redshift space distribution func- 
tion fq{r,v) (Section [L3); (3) the maximum value allowed 
for the loga rithmic derivative of the caustic amplitude A 
(Section 4.2); (4) the filling factor J^/s assumed to be constant 



over the entire interval of the caustic amplitude integration 
(Section g, equation y_3|; Section |6[ equation |2l[ ). 

Our procedure is automatic and, despite the four pa- 
rameters listed above, in essence non-parametric. In fact, q 
has been kept fixed throughout our analysis and the result 
are little affected by its variations. The maximum logarith- 
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mic derivative of A also has been kept fixed to a rather large 
value which is rarely reached. Finally, the choice J-fs = 1/2 
is a consequence of the assumption of the validity of the 
hierarchical clustering scenario. The parameter k is deter- 
mined by an automatic procedure. However, particularly un- 
fortunate situations may require a tuning of this parameter 
to locate the caustics accurately. The subjective choice of 
K is required more often in the rCDM model than in the 
ACDM model. The ACDM model yields redshift diagrams 
more similar to those of real clusters, where the density con- 
trast between galaxies falling into the cluster and unrelated 
galaxies is more evident. Therefore, our method is likely to 
be robust when applied to real clusters (Geller et al. 1999a). 

Vedel & Hartwick (1998) suggest an alternative pro- 
cedure to extract the information contained in redshift di- 
agrams. This procedure is based on a maximum likelihood 
technique and therefore has the disadvantage of being con- 
strained by the assumed model for the velocity and density 
profiles. In hierachical clustering scenarios, the infall region 
dynamics can be very complex and modelling these profiles 
might turn out to be a particularly difficult task. This pro- 
cedure is likely to be successful only when we average over 
many clusters. 

At redshift z ^ 0.2, methods based on weak gravita- 
tional lensing (e.g. Kaiser, Squires, & Broadhurst 1995; 
Seitz & Schneider 1996; Squires & Kaiser 1996; Lombardi 
& Berlin 1998) can also provide a mass estimate in the 
outer regions of clusters. These methods measure all the 
mass projected along the line of sight and suffer from sys- 
tematics due to contributions from the large scale structure 
(e.g. Bartelmann 1995; Reblinsky & Bartelmann 1999); 
moreover, these methods cannot be applied to nearby clus- 
ters. On the other hand, our approach measures the local 
meiss and, in principle, is not constrained to any redshift 
range. 

Measurements of galaxy redshifts of several nearby clus- 
ters are already currently available. Geller et al. (1999a) 
have used this method to measure the mass profile of Coma 
out to 10h~^ Mpc from the cluster centre. On smaller scales 
this profile encouragingly agrees with estimates based on X- 
ray observations. Applications to other nearby Abell clusters 
are currently underway. 
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APPENDIX A: LOCATING THE CLUSTER 
CENTRE 

Defining the centre of a real cluster is not trivial. Depending 
on the available data, we can define the centre as the posi- 
tion of the cD or the D galaxy, or as the position of the peak 
of the X-ray emission (see e.g. den Hartog & Katgert 1996). 
However, these definitions are not unique; clusters may con- 
tain more than one D galaxy or may have multiple X-ray 
peaks. 

Here, according to the definition of the cluster centre 
adopted in our A^'-body simulations, we wish to locate an 
observable cluster centre as close as possible to the mini- 
mum of the gravitational potential well of the cluster. We 
therefore locate the centre of the cluster with a two-step pro- 
cedure: (1) a hierarchical method identifies galaxies in the 
sample which are cluster members; (2) an adaptive kernel 
method estimates the cluster member density distribution 
projected onto the sky. We define the centre of the clus- 
ter as the peak of this distribution; the cluster centre in 
redshift space is the median of the cluster member veloc- 
ity distribution. Note that the main goal of the hierarchical 
method (step 1) is the identification of the cluster substruc- 
ture rather than the cluster centre. The determination of 
the cluster centre with step 2 is a natural by-product of the 
substructure analysis. 

We identify members of the galaxy cluster with a hi- 
erarchical cluster analysis (e.g. Miyamoto 1990). Cluster 
analysis classifies A'^ objects using a measure of similarity 
between any two objects. Cluster analysis produces a bi- 
nary tree with similarity decreasing from the leaves to the 
root. At any level of the hierarchy the binary tree provides 
a number of distinct groups: two members within the same 
groups have similarity larger than two members within two 
different groups. 

Different definitions of similarity have been applied to 
galaxy catalogs (Materne 1978; TuUy 1987; Gourgoulhon, 
Chamaraux & Fouque 1992). Recently, Serna & Gerbal 
(1996) suggested using the galaxy pairwise binding energy 
as a measure of similarity: 



-Ev 



'G- 



+ 



1 mirrii 



2 mi + nij 



(v. 



(Al) 



where G is the gravitational constant, nii, rrij, r^, Tj, Vi, Vj 
are the masses, positions and velocities of the two galaxies. 
The reliability of the relative binding energy as a measure 
of similarity has been questioned by Gurzadyan & Mazure 
(1998) who suggest a more physically motivated method of 
grouping galaxi es. I ndeed, when applied to observed sys- 
tems, equation (Al) has two incoveniences: (1) the masses 
rrii, ruj of the two galaxies are unknown; (2) projected infor- 
mation provides only three out of the six phase-space dimen- 
sions. By ignoring the unknown components of the position 
and velocity vectors, we obtain values of the binding energy 
Eij which are smaller than the Eij^s computed with the 
full six-dimensional information; moreover, the rank of E^ 
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of different pairs is not necessarily preserved. Despite these 
shortcomings, Serna & Gerbal (1996) tested this method 
with A'^-body simulations and obtained resuls which are gen- 
erally more satisfying than other grouping algorithms. More- 
over, comparison with our simulations shows that the centre 
located with this definition of similarity is reasonably close 
to the minimum of the potential well, as we require. 

Each galaxy is located by its vector in redshift space 
Si = {ai,Si,czi), where ai,5i are its celestial coordinates 
and czi is its radial velocity. Each galaxy pair defines the 
line of sight vector 1 = (si -I- Sj)/2 and the separation vector 
s — Si — Sj. Thus, we obtain the line-of-sight component n 
of the velocity difference and the separation projected onto 
the sky r-p 



si 

IT' 



/ 2 2n1/2 

Tp = (S — TT ) ' . 



In redshift space, equation (A.1) becomes 
mim-i 1 rriimi ■> 



Eij = 



-G- 



+ 7T- 



/ p 2 rrii + rrij 

To build the binary tree we proceed as follows: 



(A2) 



(A3) 



(i) each galaxy is a group G^; 

(ii) we compute the similarity between two groups G^, 
Gi, with the single linkage method: Efj,^ = Tnm{Eij} where 
Eij is the similarity between the member i £ Gp and the 
member j € G^; 

(iii) we replace the two groups with the largest similarity 
(smallest binding energy Eij) with a group Gk. The number 
of independent groups is decreased by one; 

(iv) the procedure is repeated from (ii) until we are left 
with only one independent group. 

If we want a catalogue of disjoint groups we need to 
choose a threshold where we cut the tree. This threshold is 
somewhat arbitrary (Materne 1978). Note that we use the 
binding energy as a similarity parameter, so we cannot use 
thresholds based on the luminosity density (TuUy 1987). 
The arbitrariness of the threshold is intrinsic to the long 
range of the force of gravity: clumps are less and less bound 
as we climb up the tree (i.e. we go towards the root) but 
they are never completely independent. 

In order to identify the final threshold which determines 
the members of the main cluster, we first compile a list of 
candidate thresholds; from this list we then identify the final 
threshold. 

We adopt the following procedure to compile the list 
of candidate thresholds. This procedure does not depend 
on the ac tual value of the binding energy computed with 
equation ( |Al[ ): therefore the value adopted for the galaxy 
mass is irrelevant. 

We walk through the tree along the main branch. At 
each step the main branch having Ud descendents (leaves) 
splits into two branches. When both branches contain a sub- 
stantial fraction of the parent descendents we label the par- 
ent of the largest branch as a threshold. In other words, 
whenever the smallest branch contains a number of de- 
scendents Us > fnd, where / is a free parameter (we use 
/ = 0.1), the level of the largest branch is a threshold. Thus, 
we obtain a list of thresholds where the main branch looses 
a substantial fraction of its descendents. 

We now need to identify the final threshold. If the clus- 




(<- binary tree root) threshold hst number 



Figure Al. Velocity dispersion a computed with the leaves of 
each threshold versus the threshold number from the candidate 
list. The root of the binary tree is on the left of the abscissa axis. 
The lines referto four different lines of sight of the ACDM cluster 
shown in Fig. h . The velocity dispersion a decreases substantially 
as we move away from the root, reaches a plateau and decreases 
once again. The starting point of the plateau (filled circles) iden- 
tifies the final threshold whose leaves are the cluster members. 



ter is not completely relaxed, galaxies within the dark mat- 
ter halo of the cluster may be distributed within subclumps. 
Ai'-body simulations (e.g. Tormen, Diaferio & Syer 1998) 
suggest that the dark matter halos of these subclumps are 
easily disrupted within the main halo, thus galaxies roughly 
feel the same potential well, regardless of their parent sub- 
clump. Therefore, we do not expect the velocity dispersion a 
to vary substantially when different subclumps are included 
or excluded from the computation of a. We expect a sub- 
stantial increase of a when we include galaxies outside the 
main halo or a decrease of a when we consider the very 
central galaxies of a given subclump. 

In fact, in our simulated galaxy clusters, a computed 
with the leaves of each threshold from our candidate list de- 
creases along the main branch (from the root to the leaves), 
reaches a plateau and decreases once again when the main 
branch splits into the very internal substructure of the clus- 
ter. The most exter nal t hreshold of the plateau defines our 
final threshold. 



Fig. 



Al 



shows some examples of this trend 
of a. 

The leaves of the final threshold are the members of the 
main cluster. Once we have identified the cluster, we can go 
down in the hierarchy of the main cluster and find substruc- 
ture easily. The level of the hierarchy where substructure is 
still present gives the degree of their relative binding energy. 

We can now use the A'^ members of the main cluster to 
locate the cluster centre. We compute the two dimensional 
density distribution projected onto the sky of the cluster 
members with the adaptive kernel method described in Sec- 
tion Ol, equations ([l5|)"([l7[). The peak of the density dis- 
tribution determines the celestial coordinates of the cluster 
centre. The median of the redshifts of the cluster members 
determines the velocity czc of the cluster. We prefer the 
median to the mean because the median is a more robust 
estimate of the central value of the velocity distribution. 
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APPENDIX B: GALAXY COORDINATES IN 
THE REDSHIFT DIAGRAM 

The angular separation of each galaxy from the cluster cen- 
tre is now trivially (see Fig. hi) 



1991, ApJ, 368, L23 
Thomas P. A., Pearce F. 



r± 



czc 



sini/j. 



(Bl) 



The line-of-sight velocity uiqs requires a more careful 
consideration. In real systems, we cannot always separate 
the peculiar velocity from the Hubble velocity reliably. The 
velocity of the galaxy B with respect to the cluster centre 
C is v^"^^ = v^^"^ — v^'^'^ = vs — vc — -ffor, where vb and 
vc are the physical velocities with respect to the observer O 
and r is the separation vector CB (Fig. hi) . We can estimate 
only the redshift of the galaxy B, czb ~ vs • fos and the 
cluster redshift, czc = vc • foe, where the hat indicates 
a versor. Thus, v^"^ • fos — czb — czc cosip — Hgr sin'y 
where we have used the fact that vc — Horoc, because 
Iv^'^'^l <C Ho\roc\- Finally, the vector r is unknown; thus, 
we define the observable line-of-sight velocity of galaxy B 



«lo 



'B,C 



■ roB + Hor sin 7 = czb — czc cos xp. 



(B2) 



This relation yields the desired quantity v^"^^ • fos as 
long as the contribution _H'o''sin7 of the Hubble velocity 
with respect to the cluster centre is negligible. This is not 
the case when, at fixed angular separation ip, we move B 
away from A (see Fig. hh; in this case both r and I7I increase 
and uios becomes a poor measure of Vg'^^ • fos- However, 
B needs to be several Mpc away from A before the Hub- 
ble contribution becomes comparable to the escape velocity 
which we expect to be of several hundreds km s~^ in a mas- 
sive cluster (see equation 111). Moreover, as we move B away 
from A both the infall (radial) velocity and the line of sight 
component of the tangential velocity decrease, so iii^s from 
these galaxies will presumably be smaller than the dios de- 
termining the caustics; thus the measure of A{r) will not be 
affected. This problem, however, becomes increasingly seri- 
ous as we increase -;/), because the relative Hubble contribu- 
tion to Vg"^^ increases and the escape velocity determining 
the caustic decreases. 

Finally, note that the procedure adopted to locate the 
cluster centre also yields, at the chosen threshold, a number 
of groups distinct from the main cluster, besides a number of 
individual galaxies, i.e. galaxies that do not belong to any 
group. In the redshift diagram we include only individual 
galaxies or galaxies which belong to groups with tt < czum, 
where n is computed with the group centre coordinates or 
the individual galax y co ordinates and the cluster centre co- 
ordinates (equation A2); we set czn^i = 4000 km s~^. The 



group centre coordinates are computed with the same pro- 
cedure used for the main cluster. 
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